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Vorticity, Shocks and Magnetic Fields in Subsonic, ICM-like Turbulence 
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ABSTRACT 

We analyze data from high resolution simulations of the generation of compressible, 
magnetohydrodynamic (MHD) turbulence with properties chosen to resemble condi¬ 
tions in galaxy clusters. In particular, the flow is driven to have turbulence Mach 
number M t ~ 1/2 in an isothermal medium with an initially very weak, uniform seed 
magnetic field (j3 = Pg/Ps = 10 6 ). Since cluster turbulence is likely to result from a 
mix of sheared (solenoidal) and compressive forcing processes, we examine the distinct 
turbulence properties for both cases. In one set of simulations velocity forcing is entirely 
solenoidal (V • Su = 0), while in the other it is entirely compressive (V x Su = 0). Both 
cases develop a mixture of solenoidal and compressive turbulent motions, since each gen¬ 
erates the other. The development of compressive turbulent motions leads to shocks, 
even when the turbulence is solenoidally forced and subsonic. Shocks, in turn, produce 
and amplify vorticity, which is especially important in compressively forced turbulence. 
To clarify those processes we include a pair of appendices that look in detail at vorticity 
evolution in association with shocks. From our simulation analyses we find that mag¬ 
netic fields amplified to near saturation levels in predominantly solenoidal turbulence 
can actually enhance vorticity on small scales by concentrating and stabilizing shear. 
The properties, evolution rates and relative contributions of the kinetic and magnetic 
turbulent elements depend strongly on the character of the forcing. Specifically, shocks 
are stronger, but vorticity evolution and magnetic field amplification are slower and 
weaker when the turbulence is compressively forced. We identify a simple relation to 
estimate characteristic shock strengths in terms of the turbulence Mach number and 
the character of the forcing. Our results will be helpful in understanding flow motions 
in galaxy clusters. 
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— turbulence 
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1. Introduction 


Most of the baryonic matter in galaxy clusters and cosmic filaments exists in the form of very 
diffuse plasma. These media, both the hotter, intracluster medium (ICM) and cooler, Warm-Hot 
Intergalactic Medium (WHIM)Q are very dynamical environments with active “weather” driven 
by ongoing accretion, substructure motions, occasional, violent merger activity and at times large 
energy inputs from starburst-driven galactic winds and very fast outflows from active galaxies 


(AGNs) (Brunetti & Jones 2014, and references therein). Simulations predict, consequently, that 
ICMs contain large-scale flows (“ICM winds”) at a fair fraction of the local sound speed, and both 
simulations and observations indicate that they contain weak-to-moderate-strength shocks, contact 
discontinuities (known in clusters as “cold fronts”) and strong bulk shear. Such flows should become 


turbulent. ICM turbulence is manifest in simulations of cosmic structure formation (e.g., 

Ryu et al. 

2003, Fujita et al. 2004| Vazza et al 2009J Ruszkowski & Oh 2011; Xu et al. 2011; Miniati 2014| Vazza 

et al 2014) and there is observational evidence in clusters that also supports this (e.g., 

Churazov et 

to 

o 

o 

Schuecker et al. 2004, 

Bonafede et al. 2010 

Churazov et al. 2012). An important feature of 


the events driving ICM turbulence is that they involve both strong shear and strong compressions 
that are not necessarily coincident. Accordingly, ICM turbulence characteristics are unlikely to be 
uniform on cluster scales. 

ICM turbulence has important consequences, including non-thermal pressure support, entropy 
and metal redistribution and probably cosmic ray acceleration. In addition, since the ICMs are 
highly ionized and conducting, such turbulence may amplify even very weak seed magnetic fields 


through the small-scale, turbulent dynamo (e.g., Schekochihin et al. 2005 Ryu et al 


2005 

Ryu et al. 

2008, 

Cho 


2014). Those magnetic fields, in turn, will determine important microphysical properties of the 


ICM, including electrical and thermal conductions and kinetic viscosity that are dependent on the 
field structure and strength (e.g., Narayan fe MedvedevpOOl Kunz et al.||2012 Parrish et al.|[2012 


Gaspari & Churazov 2013| Kunz et al. 2014; Howes 2015). The possible origins and structures of 


ICM seed fields, are varied. There are many candidate sources, ranging from primordial to plasma- 


physical, ICM- and galaxy-based (e.g., 

Widrow et al. 

2012 

Ryu et al 

2012 

for reviews). In most 

cases the distributed seed held strength available for amplihcation in 

the ICM should be several 


orders of magnitude less than the ^uGauss field values inferred to be in clusters from observations 


(e.g., Bonafede et al. 2010). Recent Planck observation of cosmic microwave background (CMB) 


anisotropies, for instance, put a constraint on the upper limit of the primordial held strength, B < 
a few nGauss at a scale of ~ 1 Mpc (Planck Collaboration 2015). Fermi observation of TeV Blazars, 
on the other hand, put a lower bound of ~ 10~ 16 Gauss again at a scale of ~ 1 Mpc for the void 
magnetic held strength (e.g., Neronov &; Vovk|2010 Chen et al.||2014 ). The specihc origin, strength 
and distribution of seed helds are beyond the scope of our work here. Once kinetic turbulence 
develops on small scales, the initial amplification of weak seed helds is exponential and rapid (e.g., 
Schekochihin et al.||2005 Cho et ah||2009 ), so memory of the initial held properties is largely lost. 


1 For simplicity, hereafter, both media will be labeled “ICMs”. 

















































































- 3 - 


While the ICM magnetic fields are very important on small scales, we note that p Gauss field 
strengths are dynamically unimportant on large scales in the ICM, since the associated, character¬ 
istic Alfven speeds, va = \Z2Pb/p 100 km/sec, while characteristic ICM sound speeds, c s ~ 10 3 
km/sec. As noted, large scale flow motions, and, indeed turbulent velocities are expected to be 
significant fractions of the sound speed. 


Turbulent motions in compressible media will generally include both solenoidal (w = Vxm/ 0) 
and compressive ((1/p)dp/dt = —V •«/ 0) components, where u is velocity, oj is vorticity and 
p is mass density. In subsonic turbulence, which should be the most prevalent in the ICM, the 
compressive motions are often ignored. However, since even solenoidal turbulence produces pressure 
fluctuations ~ pu 2 (Batchelor 1951), there will be density fluctuations if the sound speed is finite, 


thus, a compressive component to the turbulence (e.g., Kida k Orszag 1990; Porter et al. 1992 


2002 Cho k Lazarian|2002; Gaspari et al.|2014 ). When the turbulence velocities become more than 


10 % of the sound speed, these will lead to shocks (e.g., Kida k Orszag 1990). If the turbulence 
is forced by compressive motions, then that compressive component typically dominates, even if 
the turbulent velocities are subsonic (e.g., Federrath et al. 2011; Konstandin et al. 2012). In fact, 


shocks are common in the ICM with Mach numbers up to a few, so some contributions from 
compressive forcing are likely. It is important to understand these relationships when discussing 
ICM turbulence. 


Turbulent dynamo amplification of the ICM magnetic field also depends on the character of 
the turbulence, since amplification comes mostly from stretching via solenoidal motions. Because 
stretching rates in solenoidal turbulence are generally fastest on the smallest non-dissipative scales, 
efficient amplification of weak magnetic fields requires that the solenoidal motions driven on large 
scales cascade to much smaller scales before they are dissipated (call those lengths “id”)- That is, 
the hydrodynamical Reynolds number, R e ,H, of the turbulent flow, R e u ~ u e L e /udid ~ (A e /^d) 4 ^ 3 , 
must be large, where u e and L e are the fluid velocity and scale of the largest eddies, Ud is the 
solenoidal velocity on the scale id, and Kolmogorov velocity scalings, Ud ~ u e (ld/Le) 1 / 3 , were used 
to obtain the final relation. 


ICM plasmas, being very hot and diffuse, are only weakly collisional, in the sense that Coulomb 
scattering mean-free-paths, Ac, often exceed 10 kpc (e.g., Schekochihin k Cowley 2006), so are 
large compared to important dynamical scales. If Coulomb collisions were responsible for turbulent 
kinetic energy dissipation in the ICM, i.e, if id ~ Ac, then from the above relations, and assuming 
turbulence with L e ~ 100 kpc, we would expect R e ,H < 100, which is probably too small to 
allow significant turbulent amplification of ICM magnetic fields. On the other hand, the presence 
of a weak magnetic field enables a host of plasma-scale instabilities, such as the firehose and 
mirror instabilities, that very likely reduce the kinetic energy dissipation scales well below Ac (e.g., 
Schekochihin k Cowley|2006 Howes et al.|2008 Kunz et al.|2014 Howes|2015 ). Thus, it is probable 


that ICM R e ,H 100, as needed to facilitate turbulent magnetic field amplification. Direct ICM 
kinetic turbulence measurements on scales below a few kpc needed to confirm this explicitly do not 
yet exist. There are, on the other hand, recent X-ray-based measurements of density fluctuations in 
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several ICMs that are consistent with Kolmogorov, inertial turbulent velocity spectra below 10 kpc 
( |Gaspari &; Churazov 2013 Zhuravleva et al. 2015), seemingly requiring kinetic energy dissipation 
scales below that, and potentially well below that. Our simulations described below are motivated 
by the expectation that, indeed id is commonly sub-kpc, so that turbulence with R e .H 100 is 
also common. 

In this same context we mention that, of course, resistive dissipation, or magnetic diffusion, 
also plays an important role in determining the ability of turbulent flows to amplify magnetic fields. 
Obviously, if magnetic diffusion is faster than field line stretching on the smallest scales of the kinetic 
turbulence, magnetic held amplification is not effective. This comparison is generally expressed in 
MHD through the magnetic Prandtl number of the fluid, P T) m = v/rj ~ T d,m/ T d,hi where u and r] 
are the fluid kinetic viscosity and resistivity, while T^ m and represent characteristic magnetic 
and hydodynamical dissipation time scales. If both viscous and resistive dissipation in ICMs were 


determined by (infrequent) Coulomb collisions; that is, if so-called Braginskii MHD (Braginskii 
1965), applied, then P rjm 10 20 (e.g., Schekochihin & Cowley 2006). This is certainly favorable 


to turbulent magnetic held amplihcation, although such an extreme Prandtl number is unlikely. 
We noted above, for a start, that ICM viscosities are likely to be much smaller than in Braginskii 
MHD. At the same time the ICM resistivity is likely to be substantially larger than the Braginskii 
MHD model. On the latter we comment that while no direct constraints on rj values appropriate 
to ICM turbulence are available, recent analyses of cluster temperature structures strongly suggest 
that ICM thermal conduction is suppressed by at least several orders of magnitude compared to 


Braginskii MHD (Gaspari &; Churazov 2013 Zhuravleva et al. 2015). Since both thermal and 


electrical transports are controlled by the effective electron mean free path, this implies that the 
electrical conductivity (resistivity) is substantially reduced (enhanced) in ICMs over values derived 
from Coulomb scattering alone. Thus, P r ,m in ICMs is likely to be considerably smaller than 
Coulomb-scattering-based estimates, although there is no basis to suggest it is less than unity. 
The most important condition in our problem is to have P r ^ m > 1 (e.g., Schekochihin et al 2007; 
Brandenburg 2014), which seems very likely. Our simulations reported here are based on numerical, 
Euler-limit MHD, where, dissipation is negligible in resolved flows, as outlined in the next section. 
Consequently, the effective P r?m ~ 1. 

An additional and very relevant constraint on the degree of magnetic field amplihcation is the 
time available. An initially weak magnetic held embedded in solenoidal turbulence for at least a 
couple of tens of large-scale eddy turnover times can be amplified to strengths approaching energy 


equipartition with the solenoidal motions (e.g. Cho et al. 2009; Federrath et al. 2011, Beresnyak 


2012). However, the available number of such large-eddy turnovers in clusters, N e , is not generally 
very large. A simple estimate would be N e ~ fd x u e /(HL e ), where H ~ 70 km/sec/Mpc is the 
Hubble parameter, u e is the largest eddy velocity, L e is the largest eddy size and fd < 1 is an 
effective time fraction on which ICM turbulence is strongly driven. For characteristic ICM values 
u e ~ several x 100 km/sec, L e ~ 100 kpc, fd ~ 0.5, the expectation for N e is of order 10 (see, also, 
e.g., Ryu et al. 2008; Cho <fc Ryu 2009). So, while it is likely that kinetic turbulence is well-developed 
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and substantial magnetic field amplification takes place, especially on relatively smaller scales, it 
is, on the other hand, less likely that clusters and cosmic filaments should develop stationary, fully 
saturated MHD turbulence states with energy equipartition between kinetic and magnetic energies 
on large-eddy scales. 


To help explore the physics of the above issues we have carried out two series of high resolution 
simulations of driven MHD turbulence in compressible, isothermal fluids. The simulations all 
begin with a very weak (/3 = P g /Ps = 10 6 ), and, for simplicity, uniform seed magnetic field in a 
stationary medium. The initial magnetic field value corresponds roughly to 0(10) nGauss levels in 
the ICM context (see the next section), although evolution of the held past very early, exponential 
growth times, is not sensitive to this choice. Analogous simulations initiated with non-uniform seed 
magnetic fields are discussed elsewhere ( Cho Sz Yoo||2012 Emerick et al.|[2015 ). Turbulence is then 
driven towards equilibrium velocities ~ (l/2)c s , where c s is the isothermal sound speed (turbulence 
Mach number Mt ~ 1/2). For characteristic ICM sound speeds, c s ~ 10 3 km/sec, the equilibrium 
turbulent velocities are then several hundred km/sec. 


Two limiting cases for turbulence forcing are presented here. In one case velocity forcing is 
entirely solenoidal (V • Su = 0), so the turbulence is predominantly solenoidal in character as it 
evolves, and the turbulent dynamo is relatively efficient. In this case vorticity in the turbulence 
(or, as it turns out, more meaningfully, the enstrophy, e = (l/2)w 2 , is initially amplified by vortex 
stretching dynamics, as in neutral fluid turbulence. But, as the magnetic held approaches satura¬ 
tion, magnetic tension forces reduce vortex stretching, while they concentrate tangential shearing 
motions inside magnetic structures, and, thereby dominate the (positive) generation of enstrophy 
in the fully MHD turbulent hows. In the other simulation case the forcing is entirely compressive, 
so there are no applied sources of vorticity. On the other hand, vorticity is seeded and amplified 
at shocks and by Maxwell stresses, but solenoidal motions remain sub-dominant to the end of the 
simulations. Then the turbulent dynamo is suppressed. In both cases shocks develop out of the 
turbulence with strengths that are roughly predictable from standard relations describing density 
fluctuation amplitudes. While obviously idealized, the two model extremes allow us to identify 
more clearly the dependencies of the turbulence properties and the resultant magnetic held growth 
on the nature of the driving. 


Since shocks play such an important role in ICM turbulence, we touch again on the fact that 
ICMs are weakly collisional plasmas, while our simulations are based on numerical approximations 


to MHD. Internal structures of collisionless shocks depend on detailed plasma processes (e.g., Wilson 


et al 2007). Yet, over sufficiently long time intervals across the full transition they necessarily 


satisfy the same jump conditions as MHD shocks, since the jump conditions are derived from 
the basic conservation laws. Similarly, our numerical MHD shocks contain structures, which, in 
this case depend on the numerical method. In order to evaluate potential influences of internal 
shock structures in our study we carried out a Favre filtering analysis of our MHD simulations and 
compared those to analytic assessments of the roles for shocks in similar flows. We found, in fact, 
good consistency between the methods, strengthening the reliability of our results in the context 









of ICM shocks. 


- 6 - 


We present and compare these two studies here. The plan of the paper is as follows. Section 
2 outlines our numerical methods and simulation parameters. In §3 we discuss necessary physics 
of vorticity (or enstrophy) and magnetic field amplification. Section 4 summarizes results of our 
analysis, while our conclusions are summarized in §5. Appendix [A] presents an analytic discussion 
of vorticity generation across shocks in order to clarify the physical basis of our results, Appendix 
[B] looks at the same issue in terms of a Favre-hltered flow analysis of our simulation results, while 
Appendix [C] outlines a novel scheme to identify and characterize shocks in grid simulations. That 
method is used here to compute area-weighted shock Mach number probability distributions. 


2. Numerical Details 


The simulations were carried out with an isothermal TVD MHD code, updated for performance 


and parallel scaling from the code presented in 

Kim et al. 

(2001 

). The code uses constrained 

transport (Dai k Woodward 1998 

Ryu et al.|1998 

) to maintain a solenoidal magnetic field (V ■ B = 


0). It does not explicitly model viscous and resistive dissipations. The cubic simulation box had 
dimensions Lq = L x = L y = L z = 10 code units with periodic boundaries. Initially the medium 
was all at rest and had a uniform density, po = 1 , gas pressure, P g fi = 1 (so isothermal sound 
speed, c s = 1 ) and a uniform magnetic field in the x-direction with /?o = P g ,o/PB,o = 10 6 , so an 
Alfven velocity, va = \PPPb[p ~ 1-4 x 10 -3 . With these inputs the box sound crossing time is 10 
code units. Since our turbulence velocities become comparable to the sound speed, and the energy 
containing scale is about 2/3 the box size (see below), we would expect a characteristic time for 
turbulence development to be 0 ( 10 ) in these units, as it was. 


Although we emphasize that our intent has been to make these simulations as scale free as 
possible, it can be helpful to some readers to have some appropriate characteristic ICM scales in 
mind. In that spirit, if we imagine that Lq ~ 100 kpc, while c s ~ 10 3 km/sec, then the equivalent 
sound crossing time would be ~ 100 Myr; i.e., a unit of time in these simulations would then roughly 
correspond to 10 Myr. Similarly, the ICM gas pressure, P g = c 2 s p ~ 2 x 10~ n n e _3 dyne/cm 2 , 
where n e - 3 is the electron density compared to 10 ^ 3 cm -3 . The initial magnetic field with /3q = 10 6 
would correspond to Bo ~ 2 x 10 -8 Gauss, while the strongest ending RMS field values in these 
simulations (Model S2K) with (/3) ~ 0.055 would correspond to Brms ~ 5 pGauss. 


As noted in the Introduction, ICM turbulence is expected to be driven by accretion, mergers, 


vorticity generated at shocks formed during the formation of the large scale structure. Yet, defining 
the nature of the source of turbulence is not trivial. Here we adopt the common approach where 
turbulence is driven on scales comparable to the box size and the resulting turbulence developed 
on smaller scales is examined. Turbulence in our simulations was driven with a method similar to 


and substructure motions, as well as galactic winds and AGN outflows (see, e.g., Brunetti k Jones 


2014). Ryu et al. (2008), for instance, proposed a scenario in which turbulence is initiated by the 
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that of 

random field with a power spectrum, P% oc A; 6 exp(— 8k/k exp ) in the interval 1 < k/ko < 10, where 
kexp/ko = 2, with &o = 27r/Lo; <5« = 5f At was added on intervals, At = 0.01Lo/c s . The forcing 
power spectrum peaks around kd ~ (3/2)&o, or around a scale, ~ (2/3)Lo- With Lq ~ 100 
kpc, for instance, ~ 67 kpc, which is of order of the scale height of a cluster core. Forcings 
have random phases, so they are temporally uncorrelated. The amplitude of the perturbations 
was tuned to produce urms = ( u 2 ) 1 / 2 ~ 1/2 (sonic Mach number A4t = urms/cs ~ 1/2) at 
saturation. In practice the equilibrium, total turbulent kinetic-energy-equivalent velocities fell in 
the range 0.4 < urms < 0-6, depending on the nature of the driving defined below and on the 
resulting final magnetic field strength. It will be convenient below to refer times to an effective 
turbulence driving time, td = Ld/uRMS■ I n these simulations, td (4/3 )L 0 /c s 13 (in code units). 

ICM turbulence should derive from dynamics that lead to a mix of solenoidal and compressive 
driving conditions. We idealized that here by using a Helmholtz decomposition of the driving 
velocity field. It was thereby separated into solenoidal (V • Su = 0) and compressive (V x 6u = 
0) components. The fraction of the total driving kinetic energy put into solenoidal motions is 
designated below by the symbol, f s . Results are presented here for purely solenoidal, f s = 1, and 
purely compressive, f s = 0, driving. The simulations were carried out in each case over a wide range 
of grid resolutions, Ax = L$/N. For the f s = 1 case we show results from simulations carried out 
on N 3 = 1024 3 and A 3 = 2048 3 cell grids. For compressive driving, f s = 0, we include results from 
simulations on N 3 = 512 3 and N 3 = 1024 3 cell grids. 

The so-called integral scale of the velocities developed in the flows, Lp u , is a standard tool to 
help characterize turbulence properties. Specifically, 

1 / o° \ cx. 

L I,U = 2^ ( ^2 p k,udk J ^2 Pk,uk~ l dk , (1) 

\k=ko J k=ko 

where Pk )U is the power spectrum for the three-dimensional (3D) velocity field. The integral scale, 

Li )U in the fully-developed, driven turbulence will typically be comparable to, but a bit smaller 
than the driving scale, Ld, mentioned above. For the simulations reported here ~ 0.7 Ld 
(see Table [I]). A characteristic eddy turnover time at the integral scale can then be defined as 
t e j = Lj u /urms = (Ll,u/Ld)td- In our simulations t e j ~ 10 (in code units). 

As mentioned above, all of the numerical models presented here were computed using a version 
of the isothermal (P g oc p) MHD TVD code in which explicit viscosity and resistivity are absent; the 
simulations were not intended to resolve the viscous or resistive dissipation scales of the flows. They 
employ, instead, nonlinear switches that are sensitive to variations in fields at the limit of what can 
be resolved on the computational grid and apply minimal damping needed to stabilize numerical 
instabilities at Euler discontinuities including shocks, contact discontinuities and slip surfaces. For 


Stone et al. 


(1999) and 


Mac Low 


(1999). Velocity forcing 5f was drawn from a Gaussian 










this reason such simulations are sometimes called “Implicit Large Eddy Simulations” (ILESf) they 
implicity assume sub-grid structures corresponding to discontinuities (see §1 and Appendix B for 
further comments on this issue) 


Although continuously driven simulations of this kind develop a Kolmogorov-like inertial range, 
they do not have physically defined kinetic or resistive dissipation scales, such as those defined in 
Navier-Stokes flow. Thus, it is not possible to establish corresponding physical hydrodynamical and 
magnetic Reynolds numbers or, for that matter the magnetic Prandtl number. However, the simula¬ 
tions do have well-defined dissipation scales related to the mesh scale. Peak to trough variations are 
typically spread over about 4Ax, corresponding to a dissipation length, i^ ~ 4Ax = 4Lq/N. This 
grid-based dissipation scale is approximately the same for all fields. Hence, the effective Prandtl 
number is effectively P rjm ~ 1. All the presented models exhibit inertial-like power spectra when 
turbulence first develops (see Figures [3] and 11), but before magnetic stresses begin to influence 
the flows. Thus, we can construct estimates for effective hydrodynamical Reynolds numbers of the 
turbulence based on the integral scale defined in equation |lj as R e ,H = (L/,u/^d) 4//3 , where the 
4/3 index assumes for simplicity Kolmogorov velocity scaling. Reynolds numbers computed in this 
fashion range between about 240 and 1390 for the simulations presented here. 

Table [l] summarizes the basic inputs and characteristics of four models discussed in the follow¬ 
ing sections. 


3. Evolution of Vorticity and Magnetic Fields 
3.1. Vorticity and Enstrophy Generation 


Vorticity is a key measure of solenoidal turbulence, since solenoidal turbulence always includes 
circulation, and vorticity provides a measure of eddy circulation rates. It is important, however, 
to keep in mind that solenoidal turbulent energy is generally concentrated on largest-eddy scales, 
(<5u 2 oc Z 2 / 3 for Kolmogorov turbulence) while the associated vorticity generally cascades to smaller 
scales closer to dissipation scales (|w/ oc Z -1 / 3 in the above case). This difference will be important 
in our discussion of results in §4. To streamline our discussion there we briefly review some of the 
basics of its generation and amplification. The following equation governing vorticity is obtained 
from the curl of the Navier Stokes equation with magnetic (Maxwell) stresses, j x B , added; 


duo 

~dt 


— V • (ulo) + (Co ■ V)ft T 


—L-V p x V-Pt + V x 
P 2 



-I r j V^cu + V x G J TV x dd rv 


(2) 


where a^rv is the impulsive driving force, Pt = P T Pb is the sum of the gas pressure, P, and the 
magnetic pressure, Pb = (1/2)P 2 in the units we use here. T = B ■ \7B is the magnetic tension, v 


2 Also see, e.g., (Aspden et al. 


20081 for further discussion of ILES. 
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is the kinematic viscosity (assumed constant), while G = (1 /p)X7p- S, with S the standard traceless 
strain tensor (e.g., Mee & Brandenburg|2006). The box-averaged mean, ( du;/dt ) = 0, since forcing 
is uncorrelated across the box in these simulations ((V x 3d rv ) = 0). We note for clarity that the 
first two RHS terms can be combined to produce the more commonly used term, V x (u X uj). 


The first RHS term in equation ([2]) accounts for conservative advection of the vector vorticity, 
while the second represents vortex stretching. The remaining terms are vorticity source or sink 
terms. Our calculations simulate isothermal fluids, so P oc p and the baroclinic source term in 
equation ©, (Vp x VP)/p 2 , vanishes everywhere. On the other hand, the magnetic pressure will 
not generally be a barotropic function of density; that is, Pb ^ Pb(p), so the magnetic field can 
contribute a quasi-baroclinic source term even in isothermal flows. However, magnetic pressure 
variations and density variations turn out to be strongly anti-correlated in turbulent flow^J so 
their gradients are, after all, nearly anti-parallel. This term is then mostly sub-dominant in our 
simulations, even when the magnetic field is not weak. The magnetic tension term in equation Q 
represents the resistance of the field to bending and stretching, which inhibits solenoidal motions, 
including those that promote vortex stretching. The magnetic terms in equation ([2]) are initially 
small in our simulations, because the initial magnetic field is weak and uniform. But, as we shall 
see, they have roles to play once motions distort the field, and especially if the magnetic field is 
strongly amplified by dynamo action. The final, non-forcing term in equation ©> accounts for 
viscous dissipation. Our simulated fluids are nominally “ideal”, since, formally, v = 0. On the 
other hand, numerical dissipation mimics viscosity on scales close to the grid resolution, producing 
effects resembling the viscous term in equation ([2]). It allows viscous dissipation, shock formation 
and, through the viscous stress (zzV X G) term, vorticity creation in shocks. Appendix [b] exam¬ 
ines vorticity evolution in mass-weighted Favre-filtered ideal flow, where we obtain an analogous 
Reynolds-stress source term in ideal flows. 


We emphasize that shocks are expected, do develop in our simulated turbulent flows and are 
important, even though the turbulence is subsonic. This will be discussed in detail in §4.1.2 and 
§4.2.2. Turbulent motions produce those shocks. The shocks, in turn, generate vorticity. Vorticity 
changes across shocks are most simply understood by looking explicitly at shock jumps. Several 
authors have explored this issue analytically (e.g., Crocco 1937; Kevlahan & Pudritz 2009). Since 
shocks turn out to be a principal generator of vorticity in our initially irrotational, compressively 
driven turbulence simulations, and since the full physics of shock generation of vorticity is not 
necessarily intuitive, we present an outline of the relevant processes in Appendix A, extended to 
include magnetic field influences. The basic consequent hydrodynamical behaviors are contained in 
equation (A7); magnetic field augmentations are contained in equation ( |A15 ). As an additional ap¬ 
proach to understanding shock generation of vorticity, our Appendix [B] examines these interactions 
in terms of filtered flows, which largely eliminates numerical artifact issues at shock discontinuities 


J This behavior reflects the well-known dominance of slow mode oscillations over fast modes in compressible MHD 
turbulence (e.g., Kowal fe Lazarian|2010). 
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in the simulations. From that analysis we find, in fact, that the unfiltered, raw analysis of our sim¬ 
ulations gives reliable results. From this we conclude in the context of ICM shocks, which involve 
weakly collisional plasmas, that the key physics associated with the vorticity evolution at shocks 
derives from conservation of mass and momentum, so is independent of the specific dissipation 
physics. 


While equation ([ 2 ]) provides an effective description of local vector vorticity evolution, u(t,x), 
the vector field is not very useful for an understanding of the global vorticity evolution in our 
periodic-box simulations, since (cJ(t)) = 0. The vorticity magnitude, |w|, or more conveniently, the 
enstrophy, e = (1/2 )cj 2 , does provide a useful measure of global evolution. The dot product of 
equation (§ with the vorticity, oj, gives an enstrophy equation, 


de 

dt 


idr T Fstretch T F n 




(3) 


where we have combined and arranged RHS “flux” and source terms to emphasize their physical 
interpretatioi^j Specifically, 


F a dv = -V • (ue) = -(eV -u + u- Ve), 

Fstretch — ^ V)u — 2e(uj * V ) U • (h, 

Fcomp = -eV -u=^~ = -V • (ue) + u-Ve, 


p dt 


UJ 


Fbaroc = — ■ (V/9 X VP), 


P- 


UJ 


Fmag = ~2 ‘ ( V P X VPfi) + UJ ■ V X [ — ] , 


~P 


Fdiss ~ VLd • ( V LO + V X G 


(4) 


where u is the unit vector in the direction of ti. 


The enstrophy advection term, F a( i v is conservative, so that its integral over our periodic 
simulation volume always vanishes. Similarly, in our simulations the baroclinic term, Fbaroc = 0. 
The magnetic term, F may , is, as noted above, initially very small, because of the initially weak 
and uniform magnetic fields. In the solenoidal driving simulations, however, this term eventually 
plays a significant role in the evolution of the enstrophy and the turbulent energy. F mag , being 
dominated by magnetic tension, generally reduces the rate of vortex stretching. It can, however, also 
concentrate and stabilize tangential shear on small scales, so lead to a net increase in enstrophy 
above the neutral fluid level. Recall from above that in Kolmogorov turbulence e; oc P 2 / 3 , so 
that the enstrophy is concentrated on small scales even without this influence from the magnetic 
tension. The remaining contributors from equation ([ 3 ]) to enstrophy increases in our simulations 


4 We drop forcing contributions in equation (31, since, the velocity forcing, Sdrv, (which is only on large scales) is 
essentially uncorrelated with Co. The net forcing contribution to de/dt, while not exactly zero, is negligibly small. 
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are the F str et c h and F comp terms that account for vortex stretching and net enstrophy production 
in compressions, respectively. Note that the fluid compression rate, — V • it, enters into both of 
F a( j v and F comp terms. However, whereas F a dv always integrates to zero in a periodic domain, 
Fcomp does not, if there is a net alignment of the velocity and enstrophy gradient fields, so that 
j F comp dV = f u ■ Ve dV > 0. This alignment is usually true in shocks (see Appendices |A| 
and [B]) , so this term turns out in our isothermal flows to represent the dominant contribution of 
shocks to positive enstrophy growth. As we see below, the F comp term is especially important 
in the compressively forced turbulence, but less so in the solenoidally forced turbulence. The 
dissipative term, Fdi ss , should also be present through numerical viscosity, although we find in our 
simulations that the other terms alone can account reasonably well for the enstrophy evolution in 
these simulations. 


3.2. Magnetic Field Amplification 


The magnetic field evolution is governed by the induction equation, whose structure, of course, 
is very similar to equation governing vorticity. For a generalized Ohm’s law in the MHD 
approximation the induction equation is (e.g., Kulsrud, et al. 1997 Boyd &; Sanderson||2003 ), 

= V x (u x B) + rjS7 2 B - -J= —-Vn e x VP e , (5) 

ot V47ren| 


Kulsrud, et al. 1997 Boyd &; Sanderson|[ 2003 


where r/ is the resistivity (assumed constant and isotropic). The last term on the right in equation 
([5]), the so-called “Biermann battery” source term, is analogous to the baroclinic vorticity source 
term, and comes from different electron and ion mobilitie^] Here, n e and P e are the electron density 
and pressure, respectively. It could be important as a seed magnetic field generator in the early 
universe; the resulting magnetic field would be weak, for instance, (B) < 10~ 20 G in proto-clusters 
(e.g., Kulsrud, et al. 1997, Ryu et al 2012). We do not consider creation of magnetic fields in this 
work, only evolution of existing magnetic fields, so neglect that term. 


Similar to the case for vorticity, gj, equation ^ for the vector magnetic field, B, is not very 
useful for study of the global evolution of the magnetic field in our simulations. Analogous to that 
case, however, we can construct an equation for magnetic pressure (or equivalently energy density), 


dP B 

dt 


—V • ( uP B ) + 2 P b (B ■ V)u • B — P b \7 ■ u + r]B • V 2 B, 


( 6 ) 


where we have dropped the Biermann source term, and B is the unit vector in the direction 
of B. Analogous to equation the first RHS term in equation ^ represents advection of 
magnetic energy and integrates to zero over our periodic simulation boxes. The second RHS term 
is the magnetic field analogy to the vortex stretching term. It is the essential contributor to 


5 Note that \[Fk appears here, since we chose to express the magnetic field in units such that Pb = (1/2 )B 2 . 
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magnetic energy growth in the turbulent dynamo. Taken by itself in steady solenoidal turbulence 
without any back-reaction, it leads to an exponential growth in the magnetic pressure, with a rate, 
T ~ \B ■ Vu| ~ ui/l ~ oji, with l an effective scale over which magnetic flux tubes are stretched 
by velocity fluctuations. The third term is analogous to the F comp term in equation ([3]), while the 
final term represents resistive dissipation of magnetic energy. Again our simulations are assume 
to be “ideal”; i.e., formally that r/ = 0. However, once again, numerical approximations lead to 
diffusion and dissipation of magnetic fields on grid resolution scales, so mimic the effects of a finite 
r/. Since both the numerical resistivity and the numerical viscosity represent dissipation on the 
same, grid resolution scale, the effective magnetic Prandtl number, as discussed in §2, should be 
P r ,m = v/r) ~ 1- 


Our initial magnetic field, with Pb = 10 ~®P g , is too weak to have any significant dynamical 
influence at the start of these simulations. However, through the kinematic fluctuation dynamo the 
field is quickly amplified once solenoiodal turbulent motions have developed on small scales. Once 
the dynamo is underway the field energy should be amplified exponentially ( Eb oc exp(Tt)) on a 
timescale T _1 ~ t m = l m /u m ~ 1 /u m , where l m is the minimum length scale of solenoidal velocity 
fluctuations and u m is the associated velocity fluctuation. Following Beresnyak (2012) to allow 
for turbulent flux diffusion, we can estimate more quantitatively the expected exponential-phase 
growth rate of the magnetic energy to be T ~ 0.6/t m . Once the magnetic tension on lengths l m , 
Ti m ~ PB lm /l m -, is comparable to Reynolds stresses on those scales; i.e., when PB, m /l m ~ pu{l) 2 /lm 
or when Esilm) ~ Ei<(lm ), magnetic field growth on those scales will saturate as it back-reacts on 
the motions and reduces vortex and magnetic field stretching. We will see this directly in §4.1.1. 
The field can continue to be amplified by vortex stretching motions on larger scales until it saturates 
there for similar reasons. For Kolmogorov scaling the kinetic energy described by this evolution 
scales as Ek(1) oc Z 2 / 3 . Since, with this scaling, the eddy turnover time is t e = l/u(l ) oc Z 2 / 3 , the 
magnetic energy during this phase is expected to show a linear growth in time, Pb oc t. During this 
linear growth period the energy containing scale of the magnetic field should increase with time 
as Ib oc t 3 / 2 (again assuming Kolmogorov scaling). This phase continues until the scale reaches 
the turbulence driving scale, Ld, actually ~ (1/2 — 1/3 )Ld (e.g., Cho & Ryu 2009). By the time 
the growth in magnetic energy on large scales is balanced by Ohmic dissipation on small scales, 
roughly as Eb approaches the kinetic energy contained in solenoidal turbulent motions. How long 
that takes depends on the nature of the turbulence, but, saturation from an initially weak field will 
generally span many large-scale eddy turnover times. 


The expected magnetic field evolution in the compressively forced situation is different from the 
solenoidally forced situation in important respects. First, since small scale dynamo amplification 
depends on solenoidal motions, and those are much reduced in that case, the rate of magnetic 
energy growth will be much slower. Field amplification will also saturate at much lower energy 
levels, since the available energy to amplify the field is reduced, as well. 
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4. Discussion of Results 
4.1. Purely Solenoidal Forcing: f s = 1 

We look first at the evolution of the turbulence developed from purely solenoidal forcing, 
f s = 1 , starting with an overview of the various turbulent energy components and the enstrophy, 
then examine the properties of associated density fluctuations and shocks found in our simulations. 


4-1.1. Energy and Enstrophy Evolution 

Figure [l] illustrates the kinetic and magnetic energy evolution in simulations with two different 
resolutions (models S1K and S2K). The early kinetic turbulence evolution is essentially hydrody¬ 
namic, since, as discussed in §3, the initial magnetic field is too weak to modify motions significantly. 
Our setup gives a characteristic timescale on the driving scale, Ld ~ 6.7, of td = Ld/uRMS ~ 13. 
Solenoidal kinetic energy cascades to smaller scales on an eddy turnover time, t[ ~ l/ui l/ui, 
so hydrodynamical, solenoidal turbulence is well developed by t ~ td ~ 1 /uiL d hr the f s = 1 simu¬ 
lations. The kinetic turbulent energy, Er, briefly peaks after one large-eddy turnover near t ~ 15 
and then levels out by t ~ 25 at Er ~ 0.18, corresponding to urms = 0.6. That behavior obtained 
for the f s = 1 simulations we did with grid resolutions as coarse as 256 3 . The total kinetic plus 
magnetic energy remained close to Er ~ 0.18 from t ~ 25 forward in all the f s = 1 simulations. 
Thus, since the total thermal energy is fixed by the isothermal equation of state and does not 
change in our periodic volume, the total energy stored in the system was roughly constant from 
t ~ 25 onwards. The early enstrophy evolution mirrors this behavior, as we would expect and as is 
illustrated in Figure [ 2 J The rate of change in the total enstrophjj^J de/dt, peaks just after t = 10, 
while the cascade of solenoidal kinetic energy to small scales is developing, then adjusts, so that 
after t ~ 30 it remains close to zero. We will look more closely at the enstrophy evolution below. 

Figure [I] shows that the magnetic energy evolution histories are very similar in the two f s = 1 
simulations, with a slightly greater magnetic field enhancement in the higher resolution simulation 
reflecting the effectively smaller viscosity in that case. Much of this difference is actually created 
early on, but after t ~ 15, when enstrophy on small scales is fully developed and during the period 
of exponential magnetic held growth. The difference reflects the fact that the minimum eddy scale, 
l m , is larger in the lower resolution, S1K, simulation, so the initial growth rate, T oc Z m -2 / 3 , is 
reduced accordingly. The exponential magnetic energy growth is followed, as expected, by a period 
of linear energy growth in the interval, 20 < t < 80. The transition into linear growth begins 
when Er/Er ~ few % and ends several large-scale eddy turnover times (~ td ~ 13) later when 
Eb/Er ~ 30%. Figure [ 2 ] shows during this period that the amplification of enstrophy by way of 


6 The rate of total enstrophy change is obtained from the explicit difference between total enstrophy values in 
nearby, saved time steps in the simulations. 
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vortex stretching ( F stre tch ) is greatly reduced. This is a direct consequence of the inhibition by 
magnetic tension of stretching motions. Somewhat surprisingly, however, in the same time period 
there is a comparable increase in enstrophy amplification due to magnetic tension ( F mag ), so that 
the total rate of enstrophy amplification remains almost unchanged and balanced by dissipation. 
We explain that below. By the end of the S2K simulation at t = 130 the magnetic energy to 
kinetic energy ratio is Eb/Er ~ 1/2. The magnetic energy is still growing slowly, so over very 
long times it should eventually approach close to Eb/Er ~ 2/3, the value observed in simulations 
of incompressible turbulence (e.g., Ryu et al. 2008, Cho et ah 2009). 


After t ~ 40, Er drops in response to magnetic field back-reaction towards an eventual value, 
Er ~ 0.12 (0.13), in the S2K (S1K) simulation, corresponding to urms = -M-t ~ 0.5. The 
solenoidal component contains about 93% of the kinetic energy from f ~ 20 on, so about 7% is 
compressional (not shown in the figure). While sub-dominant, compressional motions, including 
shocks, are present, however, and we discuss their properties below. 


Figure [3] and Figure [4] illustrate important changes in distribution properties of kinetic and 
magnetic energies (and associated stresses) as the magnetic field back-reaction plays an increasingly 
significant effect. Figure [3] shows the kinetic and magnetic energy power spectra of the S2K simula¬ 
tion at t = 20 and t = 130. At the earlier time the turbulence is still essentially hydrodynamic with 
an approximately Kolmogorov form, Er(E) oc A; -5 / 3 , roughly in the interval 2 < k < 50, where 
k is in units of k m i n = ko = 2 tt/Lq. On all but very small scales the magnetic energy is much 
less than the kinetic energy still. By t = 130, when Eb is near saturation, magnetic energy is at 
least comparable to kinetic energy on all but the largest scales. Magnetic tension back-reaction has 
become important, and substantial kinetic energy has been removed from a wide range of scales 
below Lrf, flattening ER(k) in the process. Note, as part of this, that £%(/c) is not reduced for 
k > 200. In fact, at t = 130 the enstrophy power ( k 2 ER{k )) on these scales actually is increased 
from what it is at t = 20 and is strongly peaked around k ~ 200 (not plotted explicitly here). This 
is a reflection of the previously noted dominant enstrophy amplification after t ~ 40 by magnetic 
tension contributions as seen in Figure [2j F mag also is strongest on these scales, as it turns out. 
Associated changes in the character of the magnetic held are illustrated in Figure [4j which vol¬ 
ume renders the magnetic held intensity at the same two times. At the earlier time the strongest 
magnetic held is broadly distributed into relatively short “worm-like” filaments, generated by the 
cylindrical sheath of high rate of strain common in vortex tubes. By the later time, however, the 
strongest held structures are longer, but more important to our discussion here, the transverse 
character of the structures has transformed from round tubes into battened, filamentary ribbons, 
as evident in the image. As we will discuss in detail elsewhere (Porter et al. 2015), the ribbon 
cross sections consist of laminations separated roughly on the dissipative scale with magnetic held 
lines along the long axis conhning transverse, tangential shear layers. The magnetic tension in the 
ribbons both confines the shear layers (so enhances enstrophy on those scales) and stabilizes them 
against Kelvin-Helmholtz instabilities. 
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4-1-2. Compressions and Shock Generation 


As outlined in the introduction, even though these turbulent flows are subsonic and predomi¬ 
nantly solenoidal in character, pressure fluctuations must produce density fluctuations. The statis¬ 
tical properties of the compressions in our f s = 1 models are characterized in Figure [5j which shows 
evolution of the density probability distribution (PDF) for the S2K simulation at t = 20,80,130. 
The form of the PDF is approximately the log-normal distribution, originally suggested for super¬ 
sonic turbulence (Vazquez-Semadeni 1994). The standard deviation of the PDF decreases in this 
case from cr ~ 0.19 at t = 20, when the kinetic turbulence Mach number, Mt ~ 0.6, to a ~ 0.15 
at t = 130, as magnetic field tension begins to become significant and reduces the kinetic turbu¬ 
lence Mach number to Mt ~ 0.5. Both results are consistent with the simple scaling relation, 
a 2 = ln(l + 6 2 A1 2 ), using b = 0.3, in agreement with 6 ~ 1/3 predictions for solenoidal forcing in 
compressible, isothermal turbulence (e.g., Federrath et al. 2009 Konstandin et al. 2012 [Molina et 


al. 2012). Thus, characteristic density and gas pressure fluctuations in the S2K simulation after 
solenoidal turbulence develops are ~ 15 — 20%. 


Shocks form in this case especially by way of collisions between compression waves. The 
strength of the density fluctuations can then provide a way to estimate a characteristic, expected 
shock strength. In isothermal flows, shock jump condition is simply, 

5 l=p L = M 2 s -l, (7) 

P 

where A4 S is the shock Mach number. If, as a crude model, we set p = er, so match the density 
PDF standard deviation to a shock jump to define the characteristic shock Mach number, A4 S)C , 


M s ,c ~ \jl + y4n(l + 6 2 A4 2 ) a ~ <x 1 + ~ 1 + ^M t - (8) 

For our solenoidal simulations with M.t = 0.5 this leads to M . SjC — 1 ~ 1/12 ~ .08. Figure [6] 
shows the probability distribution function, /(A4 S ), measured at several times for shocks in the 
S1K simulation obtained using methods introduced in Appendix C. The dotted line in the figure 
corresponds to an exponential PDF, /(A4 S ) oc exp — l)/0.08). That is, the measured shock 

distribution is determined by a characteristic Mach number, Ai s = 1.08, very consistent with 
equation ([8]). The total area of the detected shocks with Mach numbers exceeding M. s c averaged 
over time is approximately 100 code units (~ L^), corresponding to a mean shock spacing, 4 ~ L^. 
Then the simulation box, or more physically, the largest driven eddy, contains, on average, one 
such shock spanning that eddy scale. 
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4.2. Compressive Forcing: f s = 0 

4-2.1. Energy and Enstrophy Evolution 


The evolution and properties of the compressively forced turbulence with f s = 0 have several 
important differences in comparison to the solenoidally forced turbulence with f s = 1. Figure [7] 
illustrates the energy evolution of the two kinetic turbulence components (solenoidal and compres¬ 
sive) along with the magnetic energy for the C1K simulation. In this case the total turbulent 
energy is almost fully accounted for over the entire simulation time by the compressive kinetic 
energy, as one might expect. An approximate total and compressive energy steady state condition, 
with Ek ~ 0.08 (Mt = 0.4) is achieved before t = 10, so within one sound crossing time of the 
simulation box. However, after t ~ 5 there is a rapid growth in solenoidal kinetic energy, and 
that component saturates at a level Ek, s /Ek,c ~ 1/15 after t ~ 20 — 25. Magnetic energy grows in 
response to the solenoidal motions, but Eb/Ekj < 0.01 even at the end of the simulation (t = 120). 

The origins of the solenoidal kinetic energy in this simulation are made clear by examination 
of Figure [8] and by noting that significant shocks from collisions of compressively driven motions 
first form around t > 5 as in Figure |9j In the interval 5 < t < 10 the compressive kinetic energy 
approaches its equilibrium level, and shocks resulting from collisions between compressive waves 
become common. Figure [8] shows the volume-integrated enstrophy growth rate, de/dt, before 
t ~ 16, along with the volume-integrated enstrophy “flux terms”, F str et c h, F comp and — F mafi [^] in 
equation ^ (see also equation Q ) . Before t « 5 there is a tiny enstrophy seeding, due primarily 
to numerical truncation. But de/dt increases abruptly around t = 5, when shocks first develop, 
growing almost three orders of magnitude before t = 10. The total growth rate is almost perfectly 
matched during that interval by F comp . By t « 10 the solenoidal kinetic energy, Ek, s , is approaching 
its equilibrium level (Figure [7]) and magnetic tension ( F mag ) and viscous dissipation (numerical, of 
course) have begun to inhibit enstrophy growth. Consequently, the enstrophy growth rate peaks 
and is then somewhat less than F comp . In any case, since the integrated F comp comes primarily from 
shocks (see §3.1, Appendix [B]) , it is then clear that enstrophy growth during the early evolution 
of the compressively forced turbulence is a by-product of shock formation. That point is shown 
explicitly in Figure [9j which shows zoomed-in images of shock structures and associated enstrophy 
source terms in a two-dimensional (2D) slice of one of the f s = 0 simulations at t = 8. 


Since the early enstrophy growth in the f s = 0 simulations is localized to shocks, we would 
expect magnetic field amplification to be similarly distributed. Indeed, that is the case. Volume 
renderings of the magnetic field distributions in the C1K simulation are shown in Figure 10 At 


t = 20 the strongest magnetic field is obviously associated with regions of recent shock passage. 
Even much later, at t = 120, when both the compressive and solenoidal turbulence motions are 
near steady-state and broadly distributed, the distribution of magnetic field is obviously clumpy 
and absent of the ribbon-like topologies noted for the solenoidal case in Figure [4j 


7 Fmag <0, so we reverse its sign in the plot. 
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The influence of shocks also shows in the turbulent energy power spectra, as shown for the 

The compressive kinetic energy 


C1K simulation at two times, t = 15 and t = 120, in Figure 11 


spectra are almost the same at the two times and resembles a Burgers form, Ex c (k) oc k~ 2 , so 
distinctly steeper than Kolmogorov, roughly over the range 5 < k < 100. The solenoidal kinetic 
energy spectrum shows modest growth during this interval over all ranges of k. It is also distinctly 
steeper than Kolmogorov, again reflecting the shock origins of these motions. 


Since the magnetic field amplification by way of the turbulent dynamo comes primarily from 
flux-tube stretching (equation ([B])), it is essentially dependent on the solenoidal turbulence compo¬ 
nent. Because that component in the f s = 0 simulations always represents less than 10% of the 
turbulent kinetic energy, it is not surprising that the magnetic field grows much more slowly and 
remains much weaker than in the f s = 1 simulations. Recall in the S1K simulation that the linear 
growth of the magnetic field approaches saturation to Eb/Ek, s ~ 1/3 around t ~ 80 — 100. In 
contrast, for the C1K, compressive forced simulation, the magnetic field is still in the linear growth 
phase with Eb/Ek, s ~ 1/8 (Eb/Ekj ~ -008) at the end of the simulation, t = 150. We extended 
a lower resolution f s = 0 simulation, ChK to t = 625, and found even on that much longer time 
interval that the magnetic energy was still very slowly increasing in time, so had not reached a 
fully saturated level. Thus, the efficiency of magnetic held energy generation from (total) turbulent 
energy input depends on the nature of the turbulence forcing. 


4-2.2. Compressions and Shock Generation 


Given the nature of forcing in the f s = 0 simulations and the resulting dominance of com¬ 
pressive kinetic energy, we should expect the density fluctuations and their associated shocks to be 
stronger than in solenoidally forced turbulence of the same Mach number, Ai t . Indeed, as Figure 


12 shows, the density PDF in our f s = 0 simulations exhibit a log-normal form, just as for our 
fs = 1 simulations, but with standard deviations, a, about twice as large. Again, applying the 
scaling relation, a 2 = ln(l + b 2 M 2 ), with J\A t = 0.4, we find b = 0.9. This is in good agreement 
with previous prediction of b ~ 1 for compressive forcing (e.g., Konstandin et ah][2012). 


As before we can compare this measure of the distribution of density fluctuations to the 
distribution of shock strengths, which is illustrated in Figure [l3j Once again the shock PDF is well 
described by an exponential form; in this case /( M. s ) oc exp (— (M. s — 1)/0.125), corresponding to 
M. sc = 1.125. The predicted value from equation ([8]) with b = 0.9 and M.t = 0.4 is Ai s ,c = 1-16, 
which, although slightly larger than measured, represents reasonable agreement, given the simplicity 
of the model. In this case the total area of detected shocks with A4 S > A4 s ,c is about 400 = 4 Lj 2 , 
so about four times greater than for the f s = 1 case. Then the mean separation between such 
shocks is l s ~ (1 /2)L^. The higher frequency of moderate strength shocks is consistent with the 
broader density distribution and the steeper energy power spectra mentioned above. 
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5. Conclusion 


Compressible turbulence in ICMs is likely to be common at levels that require consideration 
of both solenoidal and compressive motions. Both components will be present, since processes that 
drive ICM turbulence are likely to include both solenoidal and compressive forcing, and since each 
kind of turbulence can create the other. Shocks formed in the turbulence are the central link in 
that process. However, the energy partitioning of these components, as well as their energy spectra, 
depend on the character of the forces that drive them. Since ICMs are electrically conducting, 
solenoidal turbulent motions will, by the turbulent dynamo, amplify magnetic fields that thread 
them. The rate and effectiveness of that amplification depends on the strength of the solenoidal 
turbulent component, so, also on how the turbulence is driven. We found that while shocks are 
stronger, the solenoidal component contains less energy and so magnetic field amplification is less 
efficient, when turbulence is driven with compressive forcing. 


As the magnetic field becomes dynamically important, magnetic tension will inhibit solenoidal 
motions, leading to anisotropies in the turbulence (e.g., Goldreich Sz Sridhar 1995) and reductions 
in kinetic energy. Somewhat ironically, however, the magnetic tension forces can add vorticity on 
small scales by concentrating and protecting shear layers inside (~ 2D) magnetic flux ribbons. 


We proposed a simple relationship between the Mach number of the turbulence and the Mach 
numbers of shocks that they generate derivedfrom the width of the density PDF. The relation can 
be applied to both solenoidally driven or compressively driven turbulence, although the parame¬ 
terization depends on the nature of the forcing, so it is important again to evaluate the character 
of the turbulent driving forces. 
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of UNIST (1.140035.01). This work was carried out in part using computing resources at the 
University of Minnesota Supercomputing Institute. We express our gratitude to an anonymous 
referee for helpful comments that improved the presentation of the paper. 


A. Enstrophy Generation Across Shock Jumps 


Vorticity and enstrophy generation in shocked flows are fundamental to an understanding of 
compressible turbulence. To support discussions above we summarize here two simple, analytic ways 
to establish vorticity and enstrophy generation across shock jumps. We initially ignore magnetic 
field influences, but will subsequently address those. The first, mathematically more direct approach 
begins with a form of Euler’s equation, 


Du 

~Dt 


du, 

~di 


+ (u ■ V)u 


—VP, 


91 


(Al) 
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expressed explicitly in terms of vorticity using the vector identity, u • Vtt = (1/2)Vw 2 — u x Co; 
namely, 

■ 2 ' it7D ' (A2) 


- - 1 n 2 <9li 

UX LU = -VtT + -VP 9 + - 7 TT. 

2 p M at 

Near a shock with local normal direction h = q X z, defined by two orthogonal shock tangent 


directions, q and z, equation (|A2|) can be projected along q to give, 

T 


d 

u n uo z = u z u n - — I -u‘ 


1 dPg du q 


(A3) 


p ds dt 

with d/ds = q ■ V measuring variations in this direction along the shock surface. The analogous 


projection of equation (A2) in the orthogonal tangent direction, z, looks the same, with c o q and u„ 


replacing uo z and u z , while the three RHS derivatives (now with (d/ds) = z ■ V) change sign. 


Being based on ideal flow, equation (A3) does not apply inside shocks, but by using shock 


jump conditions, it can be applied to quasi-ideal flows across shocks. We emphasize that the flows 
can be inhomogeneous and unsteady, so long as the local shock jumps are well defined in space and 
time. In particular mass and momentum conservation in hydrodynamical shocks give 


1 

Un, 2 — ~^n, 1? 
^q, 2 = ^9,1? 
Uz, 2 = U Zi i, 


(A4) 


so 


and 


u\ = u\ 


r 2 — l 2 
r 2 U n,l’ 


Pg,2 ~ Pg, 1 + 


r — 1 


~ P 1 ^ 71,1 5 


(A5) 


(A6) 


where the subscripts 1,2 correspond to upstream and downstream states measured in the local 


shock rest frame, and r = 1 + p = P 2 IP\- It is straightforward using equation (A4) to show that 


5uj n = oo n p — 0J n , 1 = 0. In addition equations (A4) guarantee that d(u q p — u q ,\)/dt = 0. Using 


equations (A4), (A5) and (A6 ) we then obtain the jump in oo z as 

1 dpi 


5iO z — U} z , 2 k-h, 1 — l^z, 1 


1 


Pi ds \ 1 -t- p 



1 + p, ds 


(A7) 


where the sound speed, c s , is used in the second term inside the parenthesis, according to barotropic 
relation dP g /ds = c^dp/ds, in order to simplify the form. Alternate forms of that term, no longer 
including parenthesis and no longer assuming anything about the equation of state, would be 
[p/(u n ppi)\dPi/ds = —(p/u n p)Du q p/Dt. There is an obvious, analogous expression for 5oo q . 


Equation (A7) is equivalent to equation (4) in Kevlahan & Pudritz (2009), who present the result 
in somewhat different form. Note that we did not assume a steady shock boundary, nor a planar 
shock, nor did we enforce energy conservation across the transition. Terms including dp/ds that 
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are apparent from insertion of equation (A5) and equation (A6) into equation (|A3[) cancel, so equa¬ 


tion! A7) also incorporates compression variation, or equivalently, Mach number variation along the 


shock surface. As pointed out by Kevlahan & Pudritz (2009) this result is, then, quite general 


The first RHS term in equation (A7) represents conservation of circulation across the shock; this 


enhancement in vorticity components tangent to the shock face just reflects the reduced downstream 
area of circulation in a plane that includes the shock normal, h. The second RHS term including 
the parentheses represents baroclinic creation of vorticity across the shock when the upstream 
conditions are inhomogeneous. In an isothermal fluid, such as in our study, 1 + p, = u 2 t /c 2 , so 


that this term vanishes. The third RHS term in equation (A7) accounts for variations in the shock 


normal speed along its surface; it represents variations in the refraction of streamlines crossing the 
shock. The variations include nonuniform upstream flows in plane shocks, but also normal velocity 
variations introduced along shock surfaces due to non-planar shock geometry. In the latter context 


it is commonly associated with vorticity creation in curved shocks or intersecting shocks (Crocco 


1937). For a uniform flow into a shock with radius of curvature, R, this term has a magnitude 


6uj 


l-i m 

1 + p R 


(A8) 


An alternate derivation of equation (A7) that makes some physical contributions to the vor¬ 


ticity jump more obvious uses terms from [u ■ V)u directly from equation (Al). Again, working 
with uj z = du n /dq — du q /dn, we see that we need to evaluate jumps in two quantities, du n /dq and 
duq/dn. To avoid confusion with notation in the previous form we write here the tangential and 
normal derivatives as d/dq = V • q and d/dn = V ■ h. The first of these is easily determined from 


the first of the jump conditions (A4) as 

du nt 2 


d /I 


dq 


dq r dq 


1 Oil. 


n, 1 


u n , i dr 
r 2 dq 


(A9) 


The second can be obtained in a few more steps by projecting equation (Al) along q on each side 
of the shock (A; = 1, 2) to produce, 


du, 


q,k 


l 


dn 




1 du q,k , du qjk | 1 dPg jk | du q , k 

2 dq Uz ' k dz p dq dt 


(A10) 


Now we can, using the second and third jump conditions in equation (A4), take the difference, 
duqp/dn — duqp/dn to obtain, 


du q p du q \ 

= r— -1- 


1 


dn 


dn 


u n ,iPi 


dP. 


9,1 


dP, 


9,2 


dq dq 


Application of the shock pressure jump, equation (A6), leads finally to 
dllq.2 


dn 


_^du q p { r-ldPgp 0 r-ldunp r-lu n pdpi 

- ^ r-N I ^ ^ 


u n ,i dr 


dn u n ppi dq 


dq 


pi dq r 2 dq 


(All) 


(A12) 
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The first RHS term in equation (A12) simply accounts for compression of tangential shear across 


the shock. The remaining terms incorporate the stresses downstream caused by upstream pressure, 
density and normal velocity variations tangential to the shock, as well as tangential variations in 
the shock compression. Note that in combination with equation (A9) to create oj z , this last term 


is canceled. The others, along with the first term in equation (A9), account fully for the vorticity 


change through the shock. Indeed, the combination du n ^/dq — du q) 2 /dn from equations (A9) and 


(A12) leads once again to bio z in equation (A7). 


The enstrophy change across the shock, 6e = (1/2 )buj 2 is simply 

5e = bu ■ uj\ + —(&*;)“ = {btu z ■ uj Zj i + bu q ■ uj q ,i} + — [(bu > z ) 2 + (<5wg) 2 ] 


(A13) 


where, since bu n = 0, we have set buj ■ uj± = bu z ■ u Zt i + bu q ■ w ?1 i, and (5u) 2 = ( 5uj q ) 2 + (bu z )- 


These can be evaluated from equation (A7) and its bu q analogy. All the terms contributing to 


5e are inherently non-negative, except for the 2nd and 3rd RHS (source) terms in equation (A7) 
when they appear inside the curly bracket in equation (A13). For example, the contribution of 
the “Crocco’s” source term to the enstrophy jump equation (A13) is — [/r 2 /(l + p)](uj x V) • u n ,i, 


which can take either sign. Note that even if uj\ = 0, the last two RHS terms in equation (A13) 
will generally be finite and positive. 

On the whole, except in cases where the incident vorticity magnitude is very small and the shock 
curvature radius is very small and negative (e.g., in some shock intersections), the total change, 


be, across shocks given by equation (A13) should be positive. Thus, we expect the integrated 
quantity f u- Vedn ~ u n ,ibe > 0 across most shocks. That is significant to our analysis of volume 
integrated enstrophy evolution. It helps explain, in particular, why the volume integral of the de/dt 
contributor, F comp = — V • ue + u- Ve, (equations ([3]), Q) provides a good measure of the enstrophy 
generation by shocks in our turbulent flows (e.g., Figure [8]). Additional insights to this relationship 
come from the filtered-flow enstrophy analysis in Appendix [Bj 

The presence of magnetic fields can also contribute to and modify vorticity (enstrophy) gen¬ 


eration in shocks. Equation (A2) can be extended to MHD by including Maxwell stress terms, 


F K = - 1 (vPb-T), 

p p 


(A14) 


to the RHS. There are two ways that the Maxwell stresses enter the problem at hand. First, Fm 


explicitly provides a contribution, bujB,zi to the RHS of equation (A7), 


= ~ 


1 ( 9Pb,2 _ 

U n ,lPl V 9s 


(i + p) q s + q • T 2 — (1 + p)q ■ T\ 


(A15) 


In addition, the jump conditions expressed in equations (A4), (A5) and (A6), will be modified 


because of the anisotopic nature of T. The latter influences are complex and depend on the nature 
of the shock (fast or slow mode). A couple of simple examples, however, are sufficient to establish 
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the character of Sojb, z - First, consider the case of a plane shock interacting with a one-dimensional 
(ID) magnetic field aligned with the shock normal, B(xt)h. Then the only change from equation 
comes from the magnetic pressure, which does not in this case change across the shock, while 


the gas pressure does change according to equation (A6). The total pressure is no longer barotropic 


even in an isothermal gas, and there is an added term to equation (A7) 

1 dB 


5u b ,z = 2 /i 


Pb 


= 2ji 


u n , i 1 dB 
B ds ’ 


u n ,ipi B ds 

where Ma = u n ,i /va is the local Alfvenic Mach number of the shock. 


(A16) 


For a second example again consider a uniform magnetic field, B, now interacting with a curved 
shock having a radius of curvature, R. Then the inclination of the field to the shock will vary with 
location along the shock, as will the magnetic pressure jump even if p is constant. In addition, 


the magnetic tension forces will modify tangential velocity jumps, adding terms to equation (A5) 
of order utu n /M\. The net result of the magnetic field can be incorporate into a term with 
magnitude, 

5u B ,z ~ (A17) 


M\ R 


Since these magnetic field-induced vorticity sources scale inversely with Alfvenic Mach number 
squared compared to hydrodynamical sources, they will be relatively small in the context of ICMs 
(e.g., Ryu et al. 2008). On the other hand, these influences would become significant when magnetic 
fields are dynamically important on scales ~ R. 


B. Filtered Flows as a Probe of Enstrophy Generation in Turbulence Shocks 

The methods used in Appendix [A] to quantify the generation of vorticity and enstrophy across 
shocks are based mathematical discontinuity at shocks with well-defined states on both sides of 
the discontinuity. On the other hand, turbulent flows and the shocks they generate are by nature 
complex, and the gradient operations involved in evaluating the critical dynamical variables are 
not, in general, analytically well-defined across shocks. In addition, ICM shocks are effectively 
collisionless, so that their internal structures are both complex and unsteady. When finite dif¬ 
ference derivatives are involved in computing shock transitions these kinds of issues are similarly 
obvious. Consequently, in simulations estimates of derivative-based quantities across the numerical 
shocks might not be reliable, especially near locations of strongly curved or intersecting shocks, 
where the results of Appendix [A] suggest the most significant shock-associated vorticity generation 
probably takes place. Similar concerns exist at slip surfaces, where, again spatial derivatives are 
not necessarily properly defined. 

A simple and effective strategy to ameliorate these issues commonly employed in turbulence 
studies is to work with properly filtered flow variables that eliminate both the mathematical and 
numerical problems. In particular, derivatives of the velocity field can be expressed in terms of 
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filtered (smoothed) values together with a well-behaved sub-filter shear or Reynolds stress. Discon¬ 
tinuities formally disappear, so all the metrics are well represented within the filtered flow. This 
provides not only a mathematical and numerical framework that is well-founded, but also, as we 
show, offers useful insights into the manner in which vorticity is generated within shocks. Here we 
introduce such an analysis and apply it, in particular, to examine the shock generation of vorticity 
in flows that are initially irrotational, Cu{t, x) = 0, and to confirm our conclusion from unfiltered 
flow analyses of our simulations that the F cornp source term defined in equation Q is predominantly 
associated with shocks and represents the dominant shock-related source of enstrophy in isothermal 
flows. This finding also reinforces the previous finding that the evolution of vorticity across shocks 
is predominantly determined by basic conservation of mass and momentum across those shocks. 


We start with the unfiltered continuity and momentum equations for a compressible, ideal 
fluid (for simplicity setting B = 0 and v = 0 and ignoring the forcing term), 


! + V.(„u)=0. 


dpu 

~dT 


+ V • {puu) = —VP. 


Given any spatial normalized convolution filter 

Q{x) = / g{x - xi)Q(xi)d 3 xi 


the corresponding Favre (or mass weighted) filter (Favre 1983) is 


Q = pQ/p. 


(Bl) 

(B2) 


(B3) 


(B4) 


The resulting equations for evolution of the filtered set of fields {p, P, u) become 


dp 

~di 


+ V • (pu) = 0, 


dpu 

~dt 


+ V • {puu) = —VP — V • f, 


where t is the sub-filter scale Reynolds stress tensor, 


r = puu — puu. 


(B5) 

(B6) 


(B7) 


The resulting equation for the Favre-filtered velocity is 


du 

dt 


+ {u ■ V)u 


VP _ V • f 
p p 


(B8) 
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The curl of equation (B 8 ) provides an equation for the evolution of the vorticity of the filtered 
velocity, u = V x u, 


du -.. ~ . ~ ~ VpxVp 

——I- (u ■ V)w = (u ; • V)u — ca(V • u) H- 35 -V x 

at p z 


V-f 


(B9) 


The concerns about spatial derivatives outlined at the start of this appendix are avoided with this 
equation, since u) represents the curl of the filtered velocity, not a filtering of V x u. 


By taking the inner product of uj with equation B9 and rearranging some terms we get an 


equation for evolution of the enstrophy density ej = ^u ) 2 of the Favre-filtered flow that is very 
similar in form to equation Q for the enstrophy of the unfiltered flow; namely, 


Fadv T F s i re i c h T Fcomp T Fbaroc T Fsfs F nr ] v T Ftoti 


(BIO) 


where now 


Fadv = -V- (ue f ) = ~(e f V ■ u + u ■ Ve/), 

Fstretch — ^ V)'U — 26j((U * V)u * CU, 

Fcomp = -e/V - it = y ^ = -V • (fte/) + ft • Ve/, (Bll) 

Fbaroc = ^ ' (Vp X VP), 

P 

Fsfs ^ 

are flux and source terms analogous to those in equation Q expressed in terms of filtered quantities, 
and 

Ftot ~ Fgfrefch + F comp + Fbaroc + F s f s , (B12) 

while remembering in our isothermal simulations Fbaroc = 0. In this derivation we have neglected 
magnetic and viscous stresses, so terms F mag and Fdiss are omitted, but a new, baroclinic-like 
term has appeared, F s f s , which represents enstrophy generation through anisotropic, sub-filter 
scale Reynolds stresses. In what follows we emphasize again that uj and e/ represent solenoidal 
properties of the filtered velocity field, not a filtering of the curl of the raw velocity held. 



Since the volume integral of F a d v vanishes in our periodic domain, we focus on the three net- 
contributor source terms in equation (Bll), designating their sum as Ft 0 t (equation (B12) with 
Fbaroc = 0). For a consistency check we have verified numerically from our simulation data that, 
within about 10 %, <9e/ (x)/dt — F a d v (x) ~ Ftot(x ) over the grid, with def(x)/dt estimated explicitly 
from differences in e/(x) constructed from data in the closest saved time steps. This means that 
enstrophy accounting including only ideal terms provides a pretty good match to the data. 


As a simple, easily constructed filter, g(r), we applied a spatial convolution based on iterating 
a top-hat filter. In particular, the top-hat filter averaged a 3x3x3 block of cells to produce the 
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filtered value in the central cell. This top-hat filter was then iterated several times, which is an 
efficient way to implement a convolution with the Gaussian form, 

1 -r a /(2 <t 2 ) 

(<7\/27r) 3 

where a = 0.826 * \/N . For N > 3 this is quite accurate, and the effective convolution function is 
impressively isotropic. In our analysis we used N = 8, which is equivalent to a Gaussian convolution 
with a = 2.336Ax, where Ax is the size of one computational cell. At the price of a modest 
smoothing in the flows we thus obtained a substantial gain in ability to extract consistent solenoidal 
measures of the simulated flows in regions of sharp spatial gradients that would, otherwise, be 
difficult to characterize. 


(B13) 


Our particular interest in this Appendix is vorticity/enstrophy generation in complex shock 
structures. So, from here we focus on the application of this formalism to the early evolution 
of enstrophy in our compressively forced, f s = 0, simulations, where complex, intersecting shocks 
develop directly from the forced compressions into an initially irrotational flow. We examined in our 
simulation data both the global and local behaviors of the filtered-flow quantities in equations (BIO) 
and (Bll). With regards to the former we found, with the exception of some very early contributions 
from F s f s in the filtered flow, that the time evolution of the volume integrated filtered-flow enstrophy 
flux and source terms in equation ( |B11[ ) agreed well with the analogous unhltered-flow flux terms 
in equation Q that are shown in Figure [8j 


In reference to Figure [8] and as noted before, the net early growth of the unfiltered enstrophy 
was due almost entirely to F comp once the first shocks formed (t ~ 2) and until solenoidal motions 
became broadly distributed (t ~ 15 — 20), when F str etch became competitive. On the other hand, 
Fcomp = — eV ■ u ^ 0 only if e / 0, so some seed enstrophy is needed for that term to contribute. In 
our isothermal flows the baroclinic source term, Fb aroc , vanishes; at early times our magnetic fields 
are very weak, so | F mag | is quite small. In the earliest unfiltered flow, before any shocks developed, 
mostly numerical truncation, noise was available to seed very small amounts of vorticity. On the 
other hand, once shocks formed the necessary seeds could come from anisotropic flow stresses that 
are represented by the z^V X G term in viscous flows, or by the analogous Reynolds stress term, 
F s f s , in the ideal, filtered flow in the present discussion. Indeed, we found in the interval 2 < t < 5 
of our f s = 0 simulations that the net contributions to enstrophy growth in the filtered flow came in 
comparable amounts from F s f s and F comp . After 1 ~ 5 F s / s dropped rapidly to near zero, leaving 
Fcomp alone as the dominant contributor to def/dt before t ~ 20, just as for the unfiltered flow. 


Close examination showed, in addition, that the enstrophy production in the filtered flow 
during the early time period was highly episodic and strongly correlated spatially with shocks, 
and with shock intersections, especially. That behavior is obvious in Figure [9j taken at t = 8 
from an f s = 0 simulation and displaying enstrophy sources in a 2D slice containing a complex of 
intersecting shocks penetrating a region of irrotational flow. Figure [hi] represents a horizontal (left- 
to-right) line cut across the middle of this complex at the same time, showing the total enstrophy 
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source term, F tnt , as well as the dominant contributing terms, F s f s and F comp . The sharp negative 
spike in F to t near x = 0 (also visible in Figure [ 9 ]) comes from F stre tch (not plotted), revealing the 
generation of shear within the forward shocks that is enhanced by stretching across the rear, oblique 
shock along this line. The region to the left is previously unshocked, so irrotational. Keeping in 
mind that the quantities shown are derived from Gaussian-filtered (broadened) flow fields, it is 
then evident that as the flow enters the shock complex it is the sub-filter Reynolds stress term, 
F s f s , that creates the seed enstrophy, but that as the flow proceeds through the shock complex 
the compressive term, F comp , becomes a comparable player, while the stretching term, F str etch, can 
contribute downstream. 


C. Determining Shock Strength Distributions in Grid Simulations 


Physical shocks are generally approximated mathematically as 2D surfaces; that is, the shock 
transition becomes vanishingly thin. In grid-based numerical simulations, on the other hand, a 
shock transition is usually distributed across several grid cells, making accurate identification of 
the shock surface, as well as a proper determination of the shock normal and shock transition 
properties challenging. Several approaches to this problem have been used before (e.g., Ryu et 


al. 2003 Vazza et al 2011). Here we outline a new post-processing algorithm that can be used 


to construct and analyze shock surfaces in simulations such as ours that lead to complex shock 
distributions. We find it to be especially flexible and also quite robust in identifying and cataloging 
weak shocks. We also find that it gives results very consistent with an updated version of the method 
in Ryu et al. (2003). Our immediate goal is to compute the area weighted probability distribution 


of shock Mach numbers. Other applications of the method, including those that especially utilize 
the information it extracts about the shape of shock surfaces, are obvious. 


We assume availability of a 3D snapshot of a simulated flow, and in particular the spatial 
distribution in cells of the velocity and density fields, u and p. In brief, the procedure consists of 
partitioning the computational volume into small cubes, sized so that one can fully enclose a numer¬ 
ical shock transition, but is unlikely to enclose multiple shocks. Cubes actually containing shocks 
are identified, the location of the shock surface within the cube, its orientation and strength are 
determined. Those surfaces can then be displayed and/or their properties analyzed and tabulated. 

The initial step is a division of the computational volume into small cubes containing n c x n c x 
n c = N c cells. The cubes should be large enough to span numerical shock transitions, but we also 
want to assume at most one shock structure exists inside a given cube. In those cases where shocks 
intersect within a cube we make the simplifying assumption that the structure can be interpreted 
in terms of a single shock surface. We have found n c = 5 (so N c = 125) to be a good practical 
choice. 


Two criteria are used to determine which cubes contain a shock surface. First, there must be 
a sufficiently strong compression rate ((1 /p)dp/dt = —V • u > 0) in some fraction of the cube, and 
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second, variations in p must be primarily in one direction. The compression rate criterion includes 
two steps. We require at least a critical minimum number of cells in a cube, N s c , with V • u < 0 
and that the average divergence in those N s cells is less than a critical value, du cr u . Specifically, 
inside each cube we compute 

N s = J2h(-V-u), (Cl) 

N c 


h(q ) 


0,9 < 0l 
1,9>0J * 


(C2) 


and compare N s to N SjC . We have found N s _ c = 8 to be a good practical choice. The average 
compression rate over these cells is 


< V • u > s = 


■ u)V ■ u, 


(C3) 


which is compared to du cr it . We have found du cr it = — 1 x (V- u) max , where (V • u) rnax measures the 
maximum rarefaction rate in the simulation box, to be a reliable and effective choice in turbulent 
flows. 


The local variation in the density, p, is characterized by its gradient. A quantitative measure of 
the mean square variations of p in different directions averaged over each cube can be constructed 
in terms of the average of the tensor product of Vp with itself; namely, 

3 

Cl ^ V p V p ^cube ^ 

1 

This is a real and symmetric matrix that has positive, real eigenvalues, Ai > A 2 > A 3 > 0, with 
Ai + A 2 + A 3 =< \Vp\ 2 >cube 1 and with orthogonal, unit eigenvectors, ei,e 2 ,e 3 , related to a as 
above. The three eigenvalues correspond to the mean square variation of p in the direction of the 
corresponding eigenvector. If ei is set to be the principle direction of variation, Ai is the largest 
eigenvalue. To the extent that the faction of mean square variation along ei 


Ai + A 2 + A 3 

is close to unity, variation in the other, orthogonal directions is negligible, and the variation of p is 
primarily along the direction of e\. Our condition for variations in p is 



(C4) 


fl > fcrit, (C 6 ) 

where f cr %t is less than but close to unity. We have found fcrit = 0.8 to be an effective choice. 

Having used the above steps to establish the presence and orientation of a shock inside a 
cube, we next measure ID density and velocity profiles along the shock normal, ei, determine the 
location of the shock surface inside the cube and calculate the shock jump conditions, including, 
for instance, shock Mach number (given an equation of state - isothermal in the present case). 
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We construct the requisite density and velocity profiles by binning their cell values against the 
coordinate s = e\ ■ x, where x is the computational grid coordinate. We use grid values inside a 
cylinder centered on the cube with a cylindrical axis aligned along e\. The radius of the cylinder 
can be varied, but we find that half the size of the cube works well, since this is small enough to 
have little variation in directions transverse to e\, but large enough to include a few dozen cells 
in the averages in the s bins. The length of the cylinder along e\ extends just beyond the cubical 
volume in order to capture the variation of density and velocity near the the edge of the cube. We 
then look for a well defined shock jump; i.e., rapid changes in the density and velocity across a few 
cells with relatively constant values on either side of this transition. Finally, we assign values on 
either side of the jump as the pre- and post-shock values of p and ft, which are used to establish 
the shock jump. If no well-defined shock jump is found, the cube is rejected as a shock container. 
If the cube is accepted, we use the jump centroid to locate the shock transition, which, with e±, 
locates and orients the shock plane within the cube. 

Once we have located the shock plane we can compute its effective area inside a cube. To 
do this we first locate the N t intersection points of the shock plane with the cube edges, Ej t . In 
a given shocked cube this will range between Ni = 3 and JVj = 6. A central point in the plane 
within the cube can be generated from a linear, vector average of all these cube edge intersections. 
Those points can be sorted into a list that runs, for example, clockwise around the central point of 
the plane. Then a sequence of point pairs can be combined with the central point to construct Ni 
triangles, whose areas can be added to determine the shock area inside the cube. Appropriate shock 
properties, such as Mach number can be assigned to the area. Additional steps to be discussed 
elsewhere can be used to establish other geometric properties of the shock surface, such as its 
curvature. 
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Table 1. Model Summaries 


Model 

fs 

Grid (IV 3 ) 

tend 

PeH 

L>i, u / Ao 

te,I 

(PB,end) 

S1K 

1.0 

1024 3 

160 

556 

0.447 

8.7 

4.4 x 10 -2 

S2K 

1.0 

2048 3 

130 

1392 

0.445 

9.0 

5.5 x 10“ 2 

ChK 

0.0 

512 3 

625 

242 

0.480 

12.3 

1.2 x 10~ 3 

C1K 

0.0 

1024 3 

120 

618 

0.484 

11.5 

7.4 x 10~ 4 


Plasma /3 values at t end are {/3 end } = 1 /(PB,end)- 


Note. — All simulations were done in a periodic box of size, Lq = 
10 in code units, using an isothermal equation of state, P g oc p , mean 
density, (p) = po = 1, mean gas pressure, ( P g ) = P g $ = 1, so sound 
speed, c s = 1. Time units are, then, 0.1Lo/c s . The initial magnetic 
pressure, Pb ,o = Pg, o/A) with /3o = 10 6 , so Alfven speed va = \J2Pb/P = 
1.4 x 10 -3 . The parameter f s measures the fraction of forcing power 
in solenoidal motions. The peak driving scale, L d ss (2/3)Lo (see §2). 
The hydrodynamical integral scale, Lj tU , is defined in equation [lj while 
the eddy turnover time at L I>U , t e j = L ItU /u RM s = ( Lp u /L d )t d , where 
t d = Ld/uRMS is the “driving time”. The effective hydrodynamic Reynolds 
number, R e ,H = with £ d = 4Ax = 4Lq/-/V. 
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Fig. 1.— Evolution of kinetic, Ek , magnetic, Eb and total, Et, energies for the S1K and S2K 
simulations of isothermal, compressible, solenoidally driven (/ s = 1) MHD turbulence. 










34 - 



Fig. 2.— Enstrophy growth rate, de/dt, and the contributing source terms for the S1K simulation 
with f s = 1. Quantities are dimension t~ 3 in code time units. 
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Fig. 3.— Power spectra, E(k), of kinetic ( Ex(k )) and magnetic ( Es{k )) energies at t = 20 and 
t = 130 in the S2K simulation with f s = 1. The long-dash line segment represents a Kolmogorov 
slope. 
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Fig. 4.— Magnetic energy density, Er. volume renderings in the SKI simulation at t = 20 (Left) 
and t = 130 (Right). “Cool” is weak; “hot” is strong. Opacities are chosen to isolate stronger 
fields. 
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Fig. 5.— PDF of density in the S1K simulation at three times. Red dot-dash curves represent 
lognormal distribution fits at t = 20 and t = 130 
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Fig. 6.— PDF of shock Mach numbers in the S1K simulation. Each of the data curves (/) 30,40 
and (/) 120,140 corresponds to an average for two indicated times. The dotted line represents an 
exponential PDF in the variable J\4 — 1 with a characteristic M S c = 1.08. 
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Fig. 7.— Evolution of kinetic energy, Ek, c and Ek, s , and magnetic energy, Eb , for the C1K 
simulation of compressively driven (f s = 0) MHD turbulence. 
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Fig. 8.— Early evolution of the enstrophy growth rate, de/dt, along with contributing terms from 
equation ([3]) for the C1K simulation. 
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Fig. 9.— 2D slice isolating intersecting shocks and enstrophy generation at t = 8 in a compressive 
forcing simulation, a) V • u, delineating the pattern of shocks, b) Total enstrophy source term, 
Fiat, of Favre-filtered flow as defined in Appendix [b] c) Reynolds-stress enstrophy source term, 
F s f s , of Favre-hltered flow, d) Compressive enstrophy source term, F comp . Flows are complex, but 
generally speaking “upstream” is to the left. Warm colors represent positive values; cool colors 
represent negative values, while green is zero. The horizontal line with tick marks (separated by 


0.5 length units) corresponds to the ID cut in Figure 14 
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Fig. 10.— Magnetic energy density, Eb, volume renderings in the C1K simulation at t = 20 (Left) 
and t = 120 (Right). “Cool” is weak; “warm” is strong. Opacities are chosen to isolate stronger 
held regions. 
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Fig. 11.— Power spectra of kinetic energy, Ex,c(k) and EK, s (k), and magnetic energy, Es(k), in 
the C1K simulation at t = 15 and t = 120. 
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Fig. 12.— PDF of density in the C1K simulation at t = 120. Red dot-dash curve represents a 
lognormal distribution fit. 
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Fig. 13.— PDF of shock Mach numbers in the C1K simulation. The data curve corresponds to 
an average for times t = 40, 50, 60. The dotted line represents an exponential PDF in the variable 
Ai s — 1 with a characteristic Ai S c = 1-125. 
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Fig. 14.— Plots of the total and dominant shock-concentrated enstrophy source terms from the 
filtered-flow enstrophy equation (BIO) along a horizontal cut through the shock complex shown in 
the 2D slice in Figure [9j 









